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Abstract 


The mean field Kuramoto model describing the synchronization of a popu¬ 
lation of phase oscillators with a bimodal frequency distribution is analyzed 
(by the method of multiple scales) near regions in its phase diagram corre¬ 
sponding to synchronization to phases with a time periodic order parameter. 
The richest behavior is found near the tricritical point were the incoherent, 
stationarily synchronized, “traveling wave” and “standing wave” phases co¬ 
exist. The behavior near the tricritical point can be extrapolated to the rest 
of the phase diagram. Direct Brownian simulation of the model confirms our 
findings. 


I. INTRODUCTION 


In recent years mathematical modeling and analysis of synchronization phenomena re¬ 
ceived increased attention because of its occurrence in quite different helds, such as solid 
state physics M. biological systems [^|-|^, chemical reactions P|, etc. These phenomena 
can be modeled in terms of populations of interacting, nonlinearly coupled oscillators as hrst 
proposed by Winfree |^. While the dynamic behavior of a small number of oscillators can 
be quite rich |^, here we are concerned with synchronization as a collective phenomenon for 
large populations of interacting oscillators p. 

A simple model put forth by Kuramoto and Sakaguchi 0,111 (see also |P]), consists of a 
population of coupled phase oscillators, 0i(t), having natural frequencies cUj distributed with 
a given probability density g{uj) 


N 


9i = uji + ^i{t) + Kij sin(6'j - 9i), i = 
i=i 


( 1 ) 


Here are independent white noise processes with expected values 
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( 2 ) 


Thus each oscillator tries to run independently at its own frequency while the coupling tends 
to synchronize it to all the others. When the coupling is sufficiently weak the oscillators run 
incoherently whereas beyond a certain threshold collective synchronization appears sponta¬ 
neously. So far, several particular prescriptions for the matrix Kij have been considered. 
For instance, = K > 0 only when \i — j\ = 1, and K^j = 0 otherwise (next-neighbor 


coupling) [|T^]; K^j = K/N > Q (mean-held coupling) [^,||; hierarchical coupling |jT^; ran¬ 
dom long-range coupling |T^-pffi| or even state dependent interactions []^. In the mean-held 
case, the model (|l|)-(H) can be written in a convenient form, dehning the (complex-valued) 
order-parameter 


N 


re 


iljj _ 


= -y 

N ^ 






(3) 


where r{t) > 0 measures the phase coherence of the oscillators, and measures the 
average phase. Then eq. (|^) reads 

ei= u;i +- Oi) + iiit), i = l,2,...,W (4) 


In the limit of inhnitely many oscillators, N —>■ oo, a nonlinear integro-differential equa¬ 
tion of the Fokker-Planck type was derived for the one-oscillator probability density, 

p{6,t,uj), 


dp 


d'^p d 


(5) 


the drift-term being given by 


v{d^ t,uj) = a; -f Krsinifj — 9), 


( 6 ) 


and the order-parameter amplitude by 


p 27 r r-\-oo 

re^P = 

Jo J—oo 


e^^p(9, t, uj)g{uj) d9 du. 


(7) 


The probability density is required to be 27r-periodic as a function of 9 and normalized 
according to 


r27T 

/ p{9, t, uj)d9 = 1. 

Jo 


( 8 ) 


Mean-held models such as those described above were studied, e.g, by Strogatz and 


Mirollo [T^ in case the frequency distribution, g{u}), has rehection symmetry, g{—oj) = g{oj) 
and it is unimodal [g{oj) is non-increasing for a; > 0]. In []^, the authors showed that 
for K smaller than a certain value Kc, the incoherent equiprobability distribution, po = 
l/(27r), is linearly stable, and linearly unstable for K > K^. As D 0-I-, the incoherence 
solution is still unstable for K > Kf, [= 2/ng{0) at D = 0], but it is neutrally stable 
for K < Kp. the whole spectrum of the equation linearized about po collapses to the 
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imaginary axis. In the non/mear stability issue was addressed, and the case of a bimodal 
frequency distribution was considered [^'(cu) is even and it has maxima at a; = icuo]. In this 
case, new bifurcations appear, and bifurcating synchronized states have been asymptotically 
constructed in the neighborhood of the bifurcation values of the coupling strength. The 
nonlinear stability properties of such solutions were also studied for the explicit discrete 
example g{uj) = |[(5(a; — cuo) + + cuo)], cf. [^]. A complete bifurcation study taking 

into account the symmetry properties of g{u}) was carried out by Crawford, |^. Similar 
results were obtained by Okuda and Kuramoto in the related case of mutual entrainment 
between populations of coupled oscillators with different frequencies The main results 
concerning linear stability of incoherence with a bimodal discrete frequency distribution 
are summarized in Fig.l (cf. Fig.l, p.319 in @]). Also, in Fig.5, p. 327 of |2^ a global 
bifurcation diagram left unresolved the full behavior of the oscillatory branch starting at 
K=4D. 

The purpose of this paper is to complete the investigation started in 


201, analyzing 


in detail (asymptotically) the solution living in the neighborhood of the tricritical point 
{K/D = A,uJo/D = 1) in the parameter space {K/D,uJo/D), Fig.l. It turns out that such 
a task is far from being merely a detail, since technical difficulties are nontrivial at all, 
and results allow to complete the conjectured diagram in Fig.4 as shown in Fig.5 below. In 
Section II, a two-time analysis for the Hopf bifurcation, already developed in [^, is revisited; 
in Section III, a multiscale analysis is performed near the tricritical point, generalizing the 

H 


asymptotic analysis earlier accomplished in 


The corresponding bifurcation equations 


have been solved recasting the problem into a general formalism due to Dangelmayr and 
Knobloch [El]. Numerical results designed to con&m the previous hndings are presented in 


Section IV, and these are summarized along with the analytical results in Section V. 


II. TWO-TIME SCALE ANALYSIS FOR THE HOPF BIFURCATION 


A. Linearized problems 


Here we revisit certain results given in PO. In the Hopf analysis conducted there, 
degeneracy of an eigenvalue of multiplicity two was overlooked, as pointed out by Crawford 
m- We will recall here the relevant points of the linear and nonlinear stability analysis 
near the line K = AD in Fig. 1 where a Hopf bifurcation from incoherence arises for an 
even discrete bimodal frequency distribution g{uj). The linearized eigenvalue problem for 
this case may be obtained by inserting p = l/(27r) + exp[At] p{9,uj) in (|) and (^, and then 
ignoring terms nonlinear in p\ 


D 






r2TT 


p(9, uj) d9 = 0. 


(9) 

( 10 ) 


It can be shown that there are two eigenvalues A, which solve the equation ||T9|| : 

9 {^) 


K 

~2 


r-\-oo 


1—00 X D -\- iu 


du = 1. 


( 11 ) 
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They are explicitly given by 


when 

s{‘^) = \ + '5(‘i' + ‘i'o)l • 


( 12 ) 


(13) 


Fig. 1 is straightforwardly constructed from (|T^). Above the dashed line, 4ujo > K, and the 
eigenvalues are complex. Each complex eigenvalue is doubly degenerate due to the reflection 
symmetry of g{uj) |^. By direct substitution into (|^), it can be checked that 


je 


hi — 


Z) T A T iuj 


,-ie 


h2 — 


D + \ 


lUJ 


(14) 


are two linearly independent eigenfunctions corresponding to the same semisimple complex 
eigenvalue A 


211. They are related by the reflection symmetry uj 


- 00 . 


e 


-e. When 


A is real, these eigenfunctions are complex conjugate of each other. The eigenvalue A is no 
longer semisimple but it still has multiplicity two 


21 


B. Two-time scale analysis 

Let us now recall how to use the method of multiple scales to construct the solution 


branches which bifurcate from incoherence at K = AD, ojq > D [20|. We dehne a small 


positive parameter e: which measures the departure from the critical value = AD by 

K = Kc-\-e^K2, 0 < e-C 1. (15) 

K 2 = ±1 has to be determined later according to the direction of the bifurcating branch 
and the scaling (|I5|) will be justihed later. The probability density p{9, t, 00 ; e) will be sought 
for according to the Ansatz 0: 


= — exp.j ^^£Vj(9,i,T) + 0(£‘‘) 

T = (K - K^t = 


(16) 

( 17 ) 


The rationale behind (0) is as follows. First of all, near K = K^, small disturbances from 
incoherence decay or grow according to the values of the factor 


exp[A(iF)f] ~ exp 


- K,)t + i 


(18) 


Here A(A') is given by (0|) with K given by (|l5|) and ojq > D. Hence A(iF) ~ ±iQ-\-e^K 2 {lT 
iD/D)/A, where D = \JoJq — D"^. This explains the appearance of the two distinguished time 
scales t and r. The exponential Ansatz (0) was introduced in motivated by the failure 
of the usual expansion of p in power series of £ for the particular model considered there. 
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For that model, an algebraic Ansatz yields a vertical bifurcating branch to all orders in £. In 
other models where the unknown p is everywhere non-negative, such an exponential Ansatz 
yields an asymptotic expansion (in e) with larger domain of validity than a purely algebraic 
Ansatz 

Inserting (0) and (Q) into the governing equations (^)-(§), we obtain the hierarchy 
(3.5a)-(3.7b) of 13: 


Cai = {dt — Ddg + u}dg)ai — KcRee cxi) = 0, 



0 , 


(19) 


C a 2 


erf 


-Kcde {ailm e cxi)} , 


r2-7T ( ' 




( 20 ) 


jC (as -h (Ti(T 2 y'j = -Kcde | ^cr 2 y'j Im e ai) 


-hailm e <72 + \ - K 2 [dr(Ti + dehn. e *^(e ai)]. 




-ie /-ie' 


r 2 -K I (j 2 


as + cti(T 2 + -^ j = 0. 


( 21 ) 


Here we have defined the following scalar product 

p 27 V /* + oo 


I pZtt r-\-oo 

{a{6^uj)^ P{6^uj)) = — / / a{6^uj)P{6^uj)g{uj)dujd6. 

ZTT Jo J-OO 


( 22 ) 


The solution of the homogeneous linear equation ([T9|) is a linear combination of 
/ = 1, 2 and the complex conjugates of these terms [the pi are given by (pi])]: 


ai = 


^+(^) 


D -|- i(fl -|- cu) 


p{Q.t+ 0 ) 


+ CC + 




D + i{^l — u) 


(23) 


where = uJq — and cc denotes the complex conjugate of the preceding term (in pH 
there was A_ = 0, A+ = A. Thus two terms were missing). This value of cxi has also zero 
mean, as a function of 9. Insertion of this equation in (|^) and (pID yields 


j[2 ^2i9 


j^2 „-2ie 


C(a._ + f) = 2.- + ITTW^) 


+Ae^^^A+A_ 


-i- D + ioo 


{D + iuy + 


I + cc 
+ cc. 


(24) 
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from which 


a + d= 2(D + tuj)A^A_ 

' 2 i2D + tu;)[{D + tcoy + nY 

j^2 g2i(nt+e) 

-I-i--I- cQ 

[D + i{fl + a;)] [2D + i{Q + a;)] 

j^2 g2i(Oi-0) 

^ [D + t{n - u;)][2D + i{n - u;)] 

which has also zero mean, as required. After lengthy but rather elementary calculations to 
evaluate the right-hand side of (pTD, this equation takes on the form 

£((T 3 CT 1 CT 2 -h = Q+(r,a;)e*^™+^^ + cc + Q_(r, + cc, (26) 

6 

where only the terms that may be resonant have been kept. It is natural to look for a 
solution of the form 

^3 + aiCT 2 + ^ = + cc + + cc. (27) 

6 

We determine P± by substitution of (^) into (PBf), 

[D + z(fi ± u)]Pi - ^ (1, P±) = Q±. 

Then we can solve for P±: 


P 4 . = 


Kc{l,P±) 


+ 


Q± 


2[P-|-i(f2 ± a;)] D + i{D±uj) 


(28) 


^From (ph^ and the reflection symmetry of g{uj), we know that |Arc(l, l/[P-|-z(12±a;)]) = 1, 
so that the scalar product of 1 with (|^) produces the following non-resonance conditions: 


( 


Qa 


{ 


D i(f2 -|- cn) 

Q- 


) = o 


D -I- i{D — u) 


) = 0, 


(29) 


where we set 


/ +00 

a{ijj)g{ijj)du}. 

-00 


(30) 


The zero mean condition is also satished automatically. Some more tedious calculations lead 
hnally to two nonlinear coupled ordinary differential equations for A+(r), A_(r): 

A+ = aAj^ — (/?|A_|^ -|- 7|A+|^) A+, 


i_ = aA_ - (/3|A+|2 + 7|A_n A_, 


(31) 
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where '= (i/hr, and 


1 iD 
" ^ 4 “ ’ 


/3 = 


D 




7 = 


K2 (4D2 + ^2) > 
DK 2 ( 9^)2 + IQu;^) • 


(32) 


This result favorably agrees with that of 0 when we set A_^_ = A, A_ = 0. The needed 
stability analysis is, consequently, a little more involved than that in 
new variables 


20|. Let us dehne the 


“=i^+r^+■<'=-1-4- 


(33) 


By using (^), we obtain the following system for u and v: 

u = 2 Re a u — Re (7 + f3) — Re (7 — f3) 

h = 2 ReQ(r; — 2 Re 7 uv. 


(34) 


Clearly, u = v or u = —v correspond to traveling wave (TW) solutions, while n = 0 
corresponds to standing wave (SW) solutions. The phase portrait corresponding to a, (3 
and 7 of (0) is easily found (see Fig. 2), and the explicit solutions are (up to, possibly, a 
constant phase shift) 


A+{t) = ^0, /i = Im a - Re a 

V Re 7 Re 7 

(or R+(r) = 0 and A_{t) as A+(r) above) in case of TW solutions, and 


(35) 


A+{t) = 31_(t) = 


Re a 


Re (7 + (3) 


V = Im a 


Im (7 + (3) 
Re (7 + (3) 


Re a 


(36) 


in case of SW solutions. Notice that both SW and TW bifurcate supercritically with 
Ikswll/'^TW > 1, as indicated in Fig. 3: Re(/3 + 7 ) and Re 7 are both positive when K 2 = 1; 
whereas the square roots in (|3^) and (|^ become pure imaginary if K 2 = — 1. This indicates 
that the bifurcating branches cannot be subcritical. From the phase portrait correspond¬ 
ing to (^ID , it follows that the SWs are always globally stable, while the TWs are unstable. 
Such result was pointed out in |^, following completely different methods, while in the 
analysis was restricted to the case = v^, and thus the TWs were erroneously found to be 
stable. 


III. MULTISCALE ANALYSIS NEAR THE TRICRITICAL POINT 

Asymptotic analysis near the tricritical point, P = {K/D = A,ujo/D = 1) in Fig.l, leads 
to the introduction of a third time-scale. In fact, near such a point. 
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K = K, + + 0{e^), = cooc + + 0{e^) {K, = AD, cooc = D), (37) 


and 


K , 1 


K 2 2 , ^ 


A± = -C + - ± -,/K^- - 16u,^„ « ± -v'8D(A'2 - 41 J 2 ). 


(38) 


This shows that, besides the basic time-scale (which is denoted by t), and the slow time 
T = eH (as in ||2^), an intermediate scale, say T = et, appears. Compare 


K 2 


Y^8i7(772 ~ 4a;2) 


~ exp ± A -- -T 


(39) 


with ([T8|) above. Conseqnently, the slightly different Ansatz 


p{0, (; w; e) = — exp ^ ^ eV, («, (, T, r) + 0(e®) 


(40) 


is needed. Inserting (^) and (^OD into the governing eqnations (^-(H) leads to the hierarchy 
below, instead of (0)-(0): 


Cai = {dt - Ddg + ujdg)ai + ADdg |lm e cti)| = 0, 


(41) 


r27r 


aidO = 0 , 


c{a 2 + ^\ = -ADde {ailm e ai)} - drcri, 


(42) 


r27r / (j2 ' 




c 


0-3 + (^ 1(^2 + = -4T> de jailm e cr 2 -f- y) + ^(T 2 + y^ Im e ai) 

+UJ2 Im e“*®(e*®', cri)'| - K2dglm cxi) - - dr ^cr2 y^ , 

^ 1 ^ 0-3 + ^ 1(^2 + dO = 0, 


( 43 ) 


( Ct| (TfcT2 af'' 


-ADde ^ c^ilm e *®(e *®', erg (Ti(T 2 y ) 













3 \ 2 

*^1 1 T™ „-idije’ _ \ I ,. T™ _ I *^1 \/ 


+ ( <^3 + 0'i0'2 j ® * (e* , (Ti) + a; 2 lni e * (e* , (J2 + ~) 


+a; 2 crilm ai)' + (cr 2 + ^ ) Im ^2 + 


ext 


-id /je' 




Here 


where 


-K 2 < Im e (T 2 + + cxilm e cxi) 


-5t ( cr 2 + ^ 1 - ( <X 3 + (T 1 CT 2 + ^ 


(xr 


/ cr| 0-3 0-2 crfA 

^ f 0-4 + CTiCTs + y + ^— + — j - 0 . 


(Q;( 6 ',a;),^( 6 ',a;))' = 


^ 27 r /*+oo 


27r 


'0 ^—00 


a(0, uj)(3{e, (j)g' {(j)deduj, 


( 44 ) 


(45) 


d'ujoi^) ~ o + '^o) ~ ~ 1^0)] (1^0 — ^Oc — D). 


(46) 


The solution of the homogeneous equation (Pl) for cxi is immediately found [ff = 0 in (p3|)]: 

(47) 


Ui = ——e ^ + cc, 


D + iuj 

plus terms which decay exponentially on the fast time scale, f, and which we will systemat¬ 
ically omit. Inserting this into equation (^2]), we obtain 


( ^2 + ^ ) = 


Ax , 7 ) 2^4^ 0/1 

—e ® + CC + ^ + cc, 


D + icj 


D + iu 


(48) 


wherefrom 


(jf _ 




^ e*® + cc + 


H 2 


(£> + iu)(2D + iu) 


e2ie + cc + + cc, (49) 

B-hzu ^ ^ 


and hence a 2 - Note that the term containing B{T,t) is the solution of the homogeneous 
equation associated to C [cf. (^Tf)]. Proceeding in a similar way, we obtain 


(J 3 + (Jicr 2 




K 2 — 4a;2 


A- 


B. 


+ 


A^ 


TT 


A. 


+ 


C{T,t) 


AD{D + iuj) {D + iuy {D + iu)^ {D + D + iu 




{D -|- iujY(2D + iuj) 


e*® -b cc -1- 


+CC -b 


2 AB - AAt (;^ + 2 dT7j) 


2ie 


{D -b iuj){2D -b iuj) 

^3g3i6» 


{D -b iuj){2D -b iuj){3D + ioo) 


+ cc, (50) 
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where C = C{T,t) has a meaning similar to that of A and B. /^From this we obtain 
and hnally, from (|4^) , a^. To obtain the leading order approximation, we only need to 
determine A{T,t). Now, ([^ ) holds provided that the nonresonance condition (needed to 
remove secular terms) 


( 


D + iu) 


,P{uj,T,t)) = 0 


holds, where P(a;, T, r) denotes the coefficient of on the right-hand side of 
(PI) turns out to be the “complex Duffing equation” 

Att —t(-^2 ~ ‘iuj2)A — —ppA = 0. 

Zi 0 


(61) 


Equation 


(52) 


Such equation, however, is not sufficient to determine A, in view of the two time scales on 
which A depends. The nonresonance condition for a^, i.e. an equation like that in (|^ 
where P{u!,T,t) now denotes the coefficient of e*’’ on the right-hand side of (^if), is the 
“linearized inhomogeneous Duffing equation” 


B’ 


TT 




2 . , 2 - 


Au;2)B-^{A^B + 2\A\^B) = 

{\A\^A), 


-2Att + 


5D 


23 

2 ^ 




(53) 


where an overbar denotes taking the complex conjugate. 

Equations (|52[) and (^3]) could be analyzed directly, e.g. by extending the Kuzmak-Luke 
method (see |^, Section 4.4), to hnd the bifurcating solutions in the vicinity of the tricritical 
point and their stability. However, we can take advantage from the already existing, rather 
comprehensive theory of amplitude equations for systems invariant under the 0(2) group of 
rotations {9 ^ 9 + ip) and reflections {9 —>■ —9, uj —»• —a;) developed by Dangelmayr and 
Knobloch in Our nonlinear Fokker-Planck problem has this symmetry, therefore the 
normal form near the tricritical point (a Takens-Bogdanov bifurcation) should be the same 
one that Dangelmayr and Knobloch studied. Equations (^2]) and (|5^) in fact can be used 
to reconstruct the scaled “normal form”: 


U” - e[cAJ' + C 2 (UU' + UU') U + c^\U\^U'] - (c 4 + c^\U\^)U = 0{e^), 


df' 


(54) 


studied by Dangelmayr and Knobloch in 
T = et is the slow scale. Setting 


[cf. their equations (3.3), p. 2480]; recall that 


U = A{T, t) + eB{T, r) = A{T, eT) + eB{T, eT) 


(55) 


in (|5^) , we obtain equations for A and B which are of the same form as (^) and (|5^ . We 
can then identify the parameters termed p, z/. A, C, D in [^, and thus M = 2C + D there, 
with our quantities 


2^K2 4(U2), 2’ 5’ 5^’ 

respectively. With these identihcations. Equation 


28 


25D’ 
becomes 


and — 


38 

25D' 


(56) 
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D 


fKo 


Utt - -^{K 2 - Auj 2 )U --\UYU = e [-^Ut 


23 

2 ^ 


W^Ut - ^m^U)T) + 0 {e\ (67) 


Note that —2AttS = 0{e^). The general analysis developed in for equation can be 
used for the present case, equation (^7]) [cf. |2^, equation (3.3)]. We make the substitution 

(58) 


t/(T; e) = RiT] 


in equation (|57|), separate real and imaginary parts, and then obtain the perturbed Hamil¬ 
tonian system 


dR V 2 25D ' 


Lt — s 


fKo 


28 


V 2 25T> 


R^]L, 


(59) 


where 


L = R^<j>T 


(60) 


is the angular momentum, and 

V = V(R)^^-2(k,-4u,,)R^-A ( 61 ) 

is the potential. This system may have the following special solutions (whose stability 

properties are also pointed out here): 

{i) The trivial solution, L = 0, R = 0, which corresponds to the incoherent probability 
density, p = l/27r. Such solution is stable for i ^2 < 0 if a ;2 > 0 and for {K 2 — 4a;2) < 0 
if a;2 < 0. 

(ii) The steady-states (SS), L = 0, R = Rq = \j^D {^2 — > 0, which exists provided 

that 002 > -^ 2 / 4 . This solution is always unstable. 

{Hi) The traveling waves (TW), L = Lq = Rl^2D (u 2 — ^^ 2 ^^ > 0, R = Rq = \\J > 

0, which exist provided that 7^2 > 0 and a ;2 > 197^2/56; these solutions bifurcate from 
the trivial solution at K 2 = a ;2 = 0. When 002 = 197^^2/56, the branch of TWs merges 
with the steady-state solution branch. This solution is always unstable. 

{iv) The standing waves (SW), L = 0, R = R{T) periodic. Such solutions have been found 
explicitly in Section 5.1 of [^]. The SWs branch off the trivial solution at K 2 = 002 = 0, 
exist for 002 > llii"2/19 > 0, and terminate by merging with a homoclinic orbit of the 
steady-state {ii) on the line 002 = 117^2/19 [see equation (5.8) of |^]. This solution is 
always stable. 
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All these results are depicted in Fig. 3 below, which corresponds to Fig. 4, IV-, in the 
general classihcation (stability diagrams) reported in [^, p.266. 

In Fig. 4 below, the bifurcation diagram relevant to the present problem with a ;2 > 0 is 


given (cf. Fig.5, IV-, in |2^, p.267). 

Note that the modulated wave solutions (in the terminology of |^^), i.e. with both L 
and R periodic functions, in general with different periods, do not appear in the problem 
studied in the present paper. 

In closing, observe that, to the leading order, equation (^0]) yields 


p{e,t]uj]e) ~ — 


Rei{<P+e) 

1 + ^- h CC 

D + tLJ 


(62) 


and hence, from (^, 




R 




It follows that 


R 

e —, w 

2D' ^ 


(63) 


(64) 


which shows that, essentially, the solution U{T;e) to equation (pTl) coincides with the conju¬ 
gate of the complex order parameter [dehned by (|^)]. For this reason, in Fig. 4 the ordinate 
can be either R or r. In Fig. 5, we depicted the global bifurcation diagram which completes 
the analogous one given in |^, cf. Fig. 5 there. 


IV. NUMERICAL RESULTS 

The goal of this section is to give numerical evidence of the theoretical results obtained 
thus far. To perform this task, we have integrated the stochastic Eq. (^) by a hrst-order Euler 
method with a time step At = 0.005. In all our simulations a population of iV = 50, 000 has 
been chosen, which is large enough to neglect hnite-size effects. 

The interesting region in the space of parameters is located above the tricritical point 
{K/D = A^ujq/D = 1). To explore this region and without loss of generality, we have 
kept hxed the strength of the noise to D = 1. Then we have set cuo = 2, and we have 
swept the phase diagram by moving the coupling constant, K, thereby Ending difierent 
behavior according to the results of the previous sections. Consistently with the figures 
depicted above, we have considered only values K > A, for which the incoherent solution 
p = l/27r is unstable. For these values of K, the (partially) synchronized SW states bifurcate 
supercritically and are stable until the SW branch disappears. In this section we define the 
order parameter (0) or (|^) in such a way that r(t) G [—1,1] and that the phase does not 
experiences jumps as it increases past odd integer multiples of tt. Then the order parameter 
which we should use to compare with the results of previous Sections is |r(f)| > 0. 

Let start the discussion considering K = 5.2 . In Fig. 6 , we can see that, after a short 
transient, the order parameter |r(f)| reaches a stable state characterized by time-periodic 
oscillations of large amplitude. Clearly, this value of the coupling constant belongs to the 
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domain of the SW solution. This periodic behavior is found as soon as K becomes larger 
than 4, but near the critical point the frequency of the oscillations is very high (recall that 
ujq = 2) and their amplitude quite small. This is why we do not depict such a behavior in 
any of the hgures. Moreover, when K = 5.2, the Fourier transform of the order parameter 
exhibits a large peak at a nonzero frequency, which corresponds to a relaxation oscillation. 
This peak slowly fades out as K decreases down towards K = 4 (near the bifurcation point 
the oscillation becomes sinusoidal). 

The opposite behavior is found for larger values of K. In comparison with the last hgure 
now the amplitude of the oscillations increases while the frequency decreases in a nontrivial 
way with the coupling constant as we can see in Fig. 7 for K = 6. The system still remains 
in the domain where the standing waves are stable. 

According to the theory, the SW solution should merge with the SS solution for values 
of K large enough. Indeed, this is what we observe in Fig. 8. In this case for K = 7 
the order parameter grows exponentially fast from the initial incoherent solution to the 
time-independent partially synchronized stationary state. The conjectured global bifurca¬ 
tion diagram of Fig. 5 suggests that there may be a region where the SW and the partially 
synchronized stationary solution are both stable. In order to detect the presence of bista¬ 
bility, it is more convenient to use a deterministic numerical method to solve the nonlinear 
Fokker-Planck equation. In fact, the Monte Carlo simulation averages over realizations of 
the noise. Then different realizations may go to different stable solutions in the bistability 
region, unless we are rather careful choosing convenient initial conditions within the basin 
of attraction of one solution, and a small enough time step. Then we need an enormous 
amount of computing time for a Monte Carlo simulation to distinguish the attractor with 
smaller basin of attraction in the bistability region. Thus we have used deterministic nu¬ 
merical simulations (hnite differences) of the nonlinear Fokker-Planck equation to obtain 
the results reported below, although we have checked that costly Monte Carlo simulations 
also yield the same results in several points of the bifurcation diagram. A direct numerical 
simulation of the nonlinear Fokker-Planck equation by hnite differences shows that for suf- 
hciently large cuq the region of bistability disappears. At cuq = 1-5 we have found a narrow 
region of bistability between SW and SS solutions which is illustrated in Fig. 9. Fig. 9(a) 
shows that different initial data evolve either to the SW or to the upper SS solution for 
K = 4.95. Fig. 9(b) illustrates the abrupt transition from a SW solution to the upper SS 
solution when K changes from 4.9597 to 4.9598. When ujq is larger, cuq = 2 as in Fig. 10, 
direct simulations show a smooth transition from SW to SS. This may correspond to having 
the turning point Ki of Fig. 5 close to the end point of the SW branch. 


V. SUMMARY 

We have used the method of multiple scales to study synchronization to oscillatory 
phases in the mean-held Kuramoto model with a bimodal frequency distribution. Near 
the Hopf bifurcation points our method recovers Crawford’s results: solution branches of 
stable standing wavlles (SW) and unstable traveling waves (TW) issue supercritically from 
the incoherent (non-synchronized) state. Near the tricritical point (where a line of Hopf 
bifurcations and a line of partially-synchronized stationary states coalesce) our multiple scale 
method recovers the normal form for symmetric Takens-Bogdanov bifurcations studied by 
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Dangelmayr and Knobloch. This study allows us to establish that the bifurcating branches 
given by the local analysis of Section II end as inhnite-period bifurcation solutions. The 
unstable TW branch terminates on the SS branch, whereas the SW branch collides with 
the homoclinic loop of the SS branch in a global bifurcation of hnite amplitude. All results 
obtained in Sections II and III above agree quantitatively, as it can be shown by asymptotic 
matching (see Appendix). Furthermore, there may be an interval of parameter values where 
SW and partially-synchronized stationary solutions are both stable. Brownian and direct 
hnite-difference simulations (Section IV) conhrm these results. 
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APPENDIX 


The bifurcation diagrams in Sections II and III agree in the sense that the correspond¬ 
ing solutions match asymptotically on some overlap domain. For instance, in case of TW 
solutions, A_|_ ^ 0, A_ = 0, one obtains from (|^) 


~ Ro exp 



+ A) 

G SATs/ 


Ro exp < iT 


2Duj2 


K2 (d Rl 

(A.l) 


where Rq is a constant to be found by asymptotic matching, and tD/Q = 0(1), 

R 2 > 0 hxed, as a ;2 —0 from above. On the other hand, near the tricritical point, it is 
shown in Section III that 




(A.2) 


Let us hx a ;2 > 0 in this equation and let K 2 —>■ 0 from above. Then 



(A.3) 


and inserting the latter into equation (|A.2|) , asymptotic matching with equation ( |A-lj) yields 


Ro — 



(A.4) 


The more involved case of the SW branch can be handled in a similar way, resorting to the 
results of reference [^ |. 
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FIGURES 



Fig. 1 

FIG. 1. Linear stability diagram for the incoherent solution pQ = l/(27r) and the discrete bi- 
modal frequency distribution, g{uj) = [5{uj—ujo)+5{uj+ujq)\/2 in the parameter space {K/D,ujq/D). 
Pq is linearly stable to the left of the lines K = 4Z), u)q > D (where Hopf bifurcations take place) 
and Kl{2D) = 1 +ujq/D‘^, ujq < D (where one eigenvalue of the linearized problem becomes zero). 
To the right of these lines, the incoherent solution is linearly unstable. At the tricritical point 
K = 4Z), loq = D, two eigenvalues become simultaneously zero. The dashed line separates the 
region where eigenvalues are real (below the line) from that where they are complex conjugate 
(above the line). 




(a) (b) 

Fig. 2 

FIG. 2. Phase planes (a) (u,v), and (b) (|A_|_|, |A-|) showing the critical points corresponding 
to traveling (TW) and standing wave (SW) solutions. 
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Fig. 3 

FIG. 3. Stability diagrams {K 2 ,uj 2 ) near the tricritical point. 



Fig. 4 

FIG. 4. Bifurcation diagram (K, R) near the tricritical point for loq > D hxed. K* is the 
coupling at which a subcritical branch of stationary solutions bifurcates from incoherence. 
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Fig. 5 

FIG. 5. Global bifurcation diagram including all stationary solution branches for loq > D fixed 
as conjectured from the information on bifurcating branches available near the tricritical point. 
The location of the turning point K = Ki depends on the actual value of u)q. The exchange of 
stabilities at the turning point is postulated, not demonstrated. Numerical simulations show that 
there is a narrow region of bistability between the SW and upper SS branches for u)q = 1.5D. This 
region vanishes for u)q = 2D. 
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Time t 

FIG. 6. Time evolution of the order parameter |r(t)| for coupling strength K = 5.2, and = 1, 
ujQ = 2. Time is measured in seconds, where one second means 200 time steps. We have considered 
as initial condition the incoherent solution p = l/27r. 
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FIG. 7. Time evolution of the order parameter \r{t)\ for a larger value of the coupling constant, 
K = 6. As in the previous case there are oscillations but now their amplitude is larger as well as 
the period. 



Time t 

FIG. 8. Time evolution of the order parameter towards the stable synchronized stationary 
solution for Wo = 2, Z) = 1 and K = 7. 
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(a) 



(b) 



FIG. 9. (a) Time evolution of the order parameter in the parameter region ujq = 1.5, D = 1 
and K = 4.95 where SW and SS solutions are both stable: different initial data evolve to one of 
these solutions, (b) Details on the end of the SW solution branch and abrupt transition to the SS 
at K = 4.9598. 
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FIG. 10. Time evolution of the order parameter in the parameter region coq = 2, D = 1 for 
K = 6.03,6.04,6.05 illustrating the smooth transition between the stable SW and upper SS solution 
branches. For this value of loq there is no bistability between SW and SS solutions. 


22 





